##### full data set
setwd("/Users/nicola/Documents/PISA/Our replication/")

mydata <- read.csv("Pisa data.csv")

####These are straight forward linear regressions for comparison with the multi-level models.


##############################################################
######### BFLPE basic model: now some code to run linear regression models for each PV

five.model.summary <- matrix(data = NA, 8, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability"),c("Estimates",round(rowMeans(five.model.summary)[1:4], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[5:8], digits = 4))))

##############################################################
######### BFLPE moderated by Ability: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + mydata[,60+i]*SchoolAverageAbility + mydata[,88+i]*SchoolAverageAbility, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_Ability <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "Linear Ability*SAA", "Quad Ability*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))


###Note below is code for linear regressions for comparison with Seaton's Tables 2 and 3. We did not have time to replicate these as multi-level models but the code in included here anyway for interest.

##############################################################
######### BFLPE moderated by Highest Occupation: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISEI + SAA_ZHISEI, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_HISEI <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "HISEI", "HISEI*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Highest Education: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISCED + SAA_ZHISCED, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_HISCED <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "HISCED", "HISCED*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Home Educational Resources: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHEDRES + SAA_ZHEDRES, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_HEDRES <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "HEDRES", "HEDRES*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Cultural Possessions: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZCULTPOSS + SAA_ZCULTPOSS, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_CULTPOSS <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "CULTPOSS", "CULTPOSS*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Extrinsic motivation: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINSTMOT + SAA_ZINSTMOT, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_INSTMOT <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "INSTMOT", "INSTMOT*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Intrinsic motivation: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINTMAT + SAA_ZINTMAT, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_INTMAT <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "INTMAT", "INTMAT*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))

##############################################################
######### BFLPE moderated by Maths Self-efficacy: linear regression only

five.model.summary <- matrix(data = NA, 12, 5)

for (i in 0:4) {
	
lm.out <- lm( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZMATHEFF + SAA_ZMATHEFF, data = mydata,  weights = mydata$NewWeight) ##note column 60 of mydata is ZPV1_MATH and column 88 is Quad_ZPV1

five.model.summary[,i+1] <- c(summary(lm.out)$coefficients[,1], summary(lm.out)$coefficients[,2])
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates

}

BFLPE_MATHEFF <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "MATHEFF", "MATHEFF*SAA"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4))))


########################################
### Results tables

BFLPE
BFLPE_Ability
BFLPE_HISEI
BFLPE_HISCED
BFLPE_HEDRES
BFLPE_CULTPOSS 
BFLPE_INSTMOT 
BFLPE_INTMAT
BFLPE_MATHEFF